Ultraviolet Quantum Emitters in Hexagonal Boron Nitride from Carbon Clusters

Ultraviolet (UV) quantum emitters in hexagonal boron nitride (hBN) have generated considerable interest due to their outstanding optical response. Recent experiments have identified a carbon impurity as a possible source of UV single-photon emission. Here, on the basis of first-principles calculations, we systematically evaluate the ability of substitutional carbon defects to develop the UV color centers in hBN. Of 17 defect configurations under consideration, we particularly emphasize the carbon ring defect (6C), for which the calculated zero-phonon line agrees well the experimental 4.1 eV emission signal. We also compare the optical properties of 6C with those of other relevant defects, thereby outlining the key differences in the emission mechanism. Our findings provide new insights into the strong response of this color center to external perturbations and pave the way to a robust identification of the particular carbon substitutional defects by spectroscopic methods.


SUPPLEMENTARY NOTE 1: Calculation details
The calculations were performed based on the spin-polarized DFT within the Kohn-Sham scheme as implemented in Vienna ab initio simulation package (VASP). 1,2 A standard projector augmented wave (PAW) formalism 3,4 is applied to accurately describe the spin density of valence electrons close to nuclei. The carbon defects were embedded in a 7 × 7 bulk supercell with 196 atoms. The atoms were fully relaxed with a plane wave cutoff energy of 450 eV until the forces acting on ions were less than 0.01 eV/Å. The Brillouin-zone was sampled by the single Γ-point scheme. The screened hybrid density functional of Heyd, Scuseria, and Ernzerhof (HSE) 5 was used to optimize the structure and calculate the electronic properties.
By changing the α parameter, we modified a part of nonlocal Hartree-Fock exchange to the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) 6 with fraction α to adjust the calculated band gap. Here, α = 0.32 was used which could reproduce the experimental band gap about 6 eV. The optimized interlayer distance was 3.37 Å with DFT-D3 method of Grimme. 7 The excited states were calculated by ∆SCF method. 8 The defect formation energies E f was calculated according to the following equation, where E q d is the total energy of hBN model with defect at q charge state and E p is the total energy of hBN layer without defect. µ C is the chemical potential of carbon and can be derived from pure graphite. For N-rich condition, the chemical potential µ N = 1/2E(N 2 ), which is half of nitrogen gas molecule. For N-poor condition, the chemical potential µ B is derived from pure bulk boron and µ BN = µ B + µ N . The Fermi level ϵ Fermi represents the chemical potential of electron reservoir and it is aligned to the valence band maximum (VBM) energy of perfect hBN, ϵ p VBM . The E corr (q) is the correction term for the charged system due to the existence of electrostatic interactions with periodic condition. The charge correction terms were computed by SXDEFECTALIGN code from Freysoldt method. 9 The charge transition levels for a given defect is defined as the formation energies of two charge states are equal with setting the position of the Fermi-level in the fundamental band gap of hBN. The calculated charge transition levels for the considered defects are shown in Supplementary Figure 1.
Carbon defects may be created by kinetic processes in experiments. One may define the binding energy of carbon defects to form larger carbon clusters as , where n is the number C 2 units in the cluster, i.e., it is a n = 1 for the dimer and n = 3 for the 6C ring. The largest is the binding energy the more probable is that the defect will agglomerate if the defect species is able to diffuse. We find that 6C ring defect has the largest binding energy among the considered carbon clusters, although, in terms per C 2 unit, as shown in Supplementary Table 1.
For the excited state calculations with post-DFT method, the 6C defect was incorporated into a flake of hBN, containing 27 boron and 27 nitrogen atoms. The dangling bonds were passivated by hydrogen atoms. The calculations with the second-order approximate coupled cluster singles and doubles model (CC2) 10 and the algebraic diagrammatic construction method, ADC(2) 11 were performed with Turbomole code. 12,13 The results of time-dependent (TD) DFT and n-electron valence state perturbation theory, NEVPT2(4,4) 14 were obtained by ORCA code. 15 In all calculations, we used cc-pVDZ basis set 16   Since the experimental PL signal features a short radiative lifetime, we only focus on those arrangements where the carbon atoms are closely packed within a single honeycomb.
This structure serves as a basic and pivotal unit in the honeycomb structure and larger defects can be regarded as a combination of such. Noteworthy, the delocalization of defect orbitals should naturally decrease the excitation energy; therefore, larger defect complexes were not considered. In Supplementary Figure 2  The representations of the defect orbitals in the gap could be identified through their behaviour under inversion symmetry. Based on the wavefuntions plotted in Fig. 1 , where the first right-hand side term refers to the orbital part and the second one is for spin part (the arrows indicate the spin directions). Here, we use the antisymmetrization operator, A|xy⟩ = 1 √ 2 (|xy⟩ − |yx⟩) for the singlet wavefunctions, and the symmetrization operator, S|xy⟩ = 1 √ 2 (|xy⟩ + |yx⟩), for the triplets. nd we prove the equations according to the horizontal (σ h ) and (σ v ) reflection symmetry operators. We should note for the single Under σ h symmetry, double prime representation should be excluded because the signs are kept as described as For simplicity, we further simplify the equations by omitting the symmetry operator and spin part of the wave functions as which can be transformed to

SUPPLEMENTARY NOTE 3: The electronic termŴ Hamiltonian
We could apply these on the electronic correlation termŴ .
One could see that ⟨xx|Ŵ |xx⟩ = ⟨yx|Ŵ |yx⟩, and these two are the E ′ states. ⟨yy|Ŵ |yy⟩ = −⟨xy|Ŵ |xy⟩, these belong to A ′ states. According to the Λ and ∆ values we get, it is possible to energetically order the four states. Here, we define Λ > 0 in singlet while Λ < 0 in triplet, and ∆ < 0. For singlet ⟨yy|Ŵ |yy⟩ > ⟨xy|Ŵ |xy⟩, and the lowest adibatic potential energy surface (APES) layer should demonstrate partial A ′ 2 while for triplet the lowest one contains A ′ 1 state because ⟨yy|Ŵ |yy⟩ < ⟨xy|Ŵ |xy⟩. Based on above expression, we can have matrix form of theŴ 22,23 aŝ As shown in main text, the energy difference between the two configurations, computed by ∆SCF, are 41 meV and 308 meV for singlets and triplets, respectively. From the CC2 results, we have computed Λ and ∆ of −168.5 and −619.5 meV for the singlets and of 393 meV and 7 meV for the triplets, respectively.In Table 3

SUPPLEMENTARY NOTE 4: Derivation of the product Jahn-Teller Hamiltonian
. With the 2 × 2 Pauli matrices: we could have the linear vibronic Hamiltonian in Eq. 8 given as where the diagonal part of this expression indicates that with ±x displacement, the energy of single determinants change their energies with constructive and destructive joint vibronic coupling strength F o ± F u ;Ĥ JT is a iso-stationary function for the APES of the JT system.
Based on these data, the electron-phonon coupling parameters are calculated as,

Jahn-Teller system and the lowest polaronic state
In the double degenerate Jahn-Teller (JT) system, the vibronic ground states in each E branch can be written as therefore the four single determinants can be expressed by The first one is for the lower layer while the second one is for the upper layer. Here we just focus on the ground state therefore the first one is considered. Then the effect of double where we consider the expansion within 40 oscillator quanta (n + m) for the coefficient parameters. We could see that for both the singlet and triplet, the JT distortion is irrelevant to the A ′ 2 state, as shown in Fig. 3 of the main text that the A ′ 2 does not have JT instability. As the symmetry breaking, the irreducible representations change their character as shown in Table 5. The E ′ would partially contribute to both A ′ states while majorly to A ′ 1 considering the JT mixture, as well as in the oscillator strengths shown in Table 6.

SUPPLEMENTARY NOTE 6: Intersystem crossing and
the non-radiative decay from |A ′ 1 ⟩ to |A ′ 2 ⟩ Beside the radiative decay, we have also explored a possibility of the non-radiative transition to the triplet manifold through the intersystem crossing (ISC). This process is mainly governed by the spin-orbit coupling (SOC), and the possible pathways are depicted in Fig.

4(d) in main text. The SOC interaction can be expressed aŝ
where λ x,y are the non-axial components, while λ z is an axial component. In particular, λ x,y couples the triplet states with the non-zero spin projections (m s = ±1) with singlets of  have the same electronic configuration |e ′′ e ′′ ⟩ only the axial part is non-vanishing. The SOC splits 3 E ′ into m s = ±1 sub-states A 1,2 and E 1,2 with the m s = 0 E x,y state. In addition, 1 A ′ 1 could also decay to 3 E ′ due to the mixture with 1 E ′ . The ISC rate from singlet to triplet can be calculated by 26 where F is the spectral function of vibrational overlap between the singlets and triplet states, and ∆E is the energy splitting between singlet and triplet levels. From the TDDFT calculations, we found the largest value of λ z of only 1.5 GHz. Given a considerable energy gap between the states, this translates into the enormously large τ ISC , and disables the ISC process in the zeroth-order (see Supplementary Figure 3). Yet, we note that the triplet manifold may be populated via a nongeminate recombination of hot charge carriers achieved by a two-photon absorption process. Another important issue is the non-radative decay from the JT distorted |A ′ 1 ⟩ state to the low lying dark |A ′ 2 ⟩ state. The defect would be a dim center if the non-radiative decay is faster than the radiative process and would bleach itself.
The non-radiative transition rate (Γ nrad ) between the two states can be calculated 27 as Here, W if is the electronic term and X if (T ) is temperature dependent phonon term. g is the equivalent energy-degenerate atomic configurations and p in is the thermal population. i and f correspond the initial and final state and here the C 2v and D 3h are used to generate the configuration coordinates. The phonon matrix |χ i |Q − Q 0 |χ f | sums up the harmonic oscillator wavefunctions that enters the non-radiative recombination process. ∆E if is energy difference between the two states and ψ is the single particle wavefunction from DFT. The electron-phonon coupling matrix can be mimic by the single particle wave function overlap between |e ′′ ox ⟩ and |e ′′ oy ⟩ together with |e ′′ ux ⟩ and |e ′′ uy ⟩, as shown in Supplementary Figure 4. The total transition rate from |A ′ 1 ⟩ is The averaged rate is 319 MHz. Also based on Γ rad Γ , where Γ is the total transition rate, we can estimate the quantum yield is about 62.7%, however, this is the optimal value without considering the temperature effect.

SUPPLEMENTARY NOTE 7: The temperature dependent quantum yield
Temperature influences the thermal distribution between the polaronic states |Ã ′ 1 ⟩ and |Ẽ ′ ⟩ and further influence the brightness, i.e., PL intensity. This process τ 23 is in picosecond range (0.17 ps at 100 K for singlet) and much faster than the non-radiative decay to low-lying |A ′ 2 ⟩ according to ℏω E e −δ/k B T . 28,29 At room temperature, 76% electrons occupy the bright |Ẽ ′ ⟩ state and this accounts for the radiative decay process. Here we evaluate the temperature effect according following equation, where δ is given in the main text. The k B is the Boltzmann constant. For the non-radiative part, the temperature effect is plotted in Supplementary Figure 4c. We assume the nonradiative decay rate remains the same for the two polaronic states Γ 31 ≈ Γ 21 . The temperature dependent quantum yield is plotted in Supplementary Figure 5. At low temperature region the brightness is close to zero due to small radiative decay and it is about 52.1% at 300 K. There is non-radiative decay from |A ′ 2 ⟩ to ground state, however, this is slow due to the large band gap. Therefore the |A ′ 2 ⟩ could be a relative long-live excited metastable state and can be further excited to conduction band by two-photon excitation.

SUPPLEMENTARY NOTE 8: Phonon properties related to carbon defects
The calculated phonon densities of states for the perfect hBN and one with the ring 6C defect are shown in Supplementary Figure 6. The defect introduces new phonon modes with the most prominent one locating at 200 meV. This peak also has largest contribution for the partial HR factor. We also find two degenerate modes at low energy range (3 meV) which are the collective motion of atoms, as shown in (g). They represent a mutual displacement of the hBN layers and are naturally missing for the monolayer configuration. The two modes are related to the D 3h symmetry while they are missing at single layer model and C N C B and 4C pair . Without these two modes, the HR factor can reduce to 1.7. This phenomena motivated us to calculate the single layer case. As expected, we get S = 1.8 with the same parameters in bulk. However, for single layer model, the GW band gap is about 7.1 eV, 27 therefore, to reproduce it the mixing parameter α should be further increased to 0.5. Under this condition, the calculated ZPL for 6C ring defect is 5.02 eV and S = 2.14.
In particular, to confirm the involvement of carbon in a colour centre, it is commonplace to use the isotopic purification method, that incorporates 13 C into the lattice of the material during its growth. Here, we determine the isotopic shift in the emission energy and sideband for the 6C defect associated by replacing 12 C with 13 C isotope. First, we calculated the sideband with 100% of 13 C isotopes and found that the HR factor reduces to 1.78, see

and 6C
To evaluate the strain effect, we apply -2% to +2% strain along zigzag and armchair directions as shown in Supplementary Figure 8 which could be mapped out in experiment. 30 A 5 × 3 √ 3 square supercell is used with C N C B and 6C considered and the calculation parameters are the same as mentioned in the Supplementary Note 1. With applied strain, the symmetry of carbon cluster (6C) reduces to C 2v so the pJT solution is not included.
For simplicity, here we just calculate the vertical excitation. For C N C B the strain is applied parallel (armchair) or perpendicular (zigzag) to the C-C dimer bond and the transition is between b 2 defect state so the excitation energy evolution demonstrates nearly linear dependence and decreases as increase of strain. Strain along zigzag direction shows smaller influence since it does not alter the C-C bond length while the strain along armchair does.
In 6C, the situation is much more complex. Generally, strain decreases the excitation energy however it is not linear response. We speculate the strain might be an effective way to distinguish the type of carbon defects, at least for C N C B and 6C. Also, we have to admit the strain we consider is the simplest case while the shear strain is not included. The  transition dipole moment evolution can be deduced from the imaginary part of dielectric function. For C N C B the dipole moment does not change much. As increase of strain on 6C, the intensity of dipole moment increases along zigzag direction while decreases along armchair direction. The lifetime would change more dramatically for strain along armchair since both the transition energy and transition dipole moment decrease.
In particular, even more spectacular effect is expected for the intensity of the optical signal. To this end, we considered the effect of the applied strain on the absorption intensity by TD-PBE0, see Supplementary Figure 9. The calculations show that for C N C B , the uniaxial strain does not alter the response density. As a result, no significant evolution of the intensity is observed. By contrast, in a D 3h configuration, the efficient absorption in the A 1 state of 6C is hindered by the presence of C 3 symmetry axes. Upon the effect of strain along the zigzag direction, this symmetry constraint is lifted. The response density primary localises along the armchair direction, while the change in the resulting transition dipole moment from 0 to 3% of strain increases the absorption intensity by more than an order of magnitude.